home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.3 / xah / aoi_hist.c < prev    next >
C/C++ Source or Header  |  1999-09-11  |  5KB  |  260 lines
  1. /* 
  2.  * aoi_hist.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * AOI_HIST(ogram).C
  11.  *
  12.  * routines to evaluate AOI histogram and related quantities
  13.  *
  14.  */
  15. #include "xah.h"
  16.  
  17. #define NBIN_LIM        50
  18. #define N_BINS            51
  19. #define BIN_WIDTH        5
  20.  
  21. #define MAX_AREA        125           /* refers to special appl */
  22. #define MIN_AREA        0
  23.  
  24.  
  25. #undef HIST_DEBUG
  26.  
  27. /*
  28.  * determine mean area
  29.  */
  30. double
  31. eval_mean (int nd, unsigned int *area)
  32. {
  33.   int i;
  34.   double mean = 0.0;
  35.  
  36.   for (i = 0; i < nd; i++) {
  37.     mean += (double) (*(area + i));
  38.   }
  39.   return (mean / (double) nd);
  40. }
  41.  
  42.  
  43. /*
  44.  * determine standard deviation of input data set
  45.  */
  46. double
  47. eval_sdev (double mean, int nd, unsigned int *area)
  48. {
  49.   int i;
  50.   double dai, var = 0.0;
  51.  
  52.   for (i = 0; i < nd; i++) {
  53.     dai = (double) (*(area + i)) - mean;
  54.     var += (dai * dai);
  55.   }
  56.   return (sqrt (var /= (double) (nd - 1)));
  57. }
  58.  
  59.  
  60. /*
  61.  * determine skewness of input data set
  62.  */
  63. double
  64. eval_skew (double mean, double sdev, int nd, unsigned int *area)
  65. {
  66.   int i;
  67.   double ai, dai, sk;
  68.  
  69.   sk = 0.0;
  70.   for (i = 0; i < nd; i++) {
  71.     ai = (double) (*(area + i));
  72.     dai = (ai - mean) / sdev;
  73.     sk += (dai * dai * dai);
  74.   }
  75.   return ((sk /= (double) nd));
  76. }
  77.  
  78.  
  79.  
  80. /*
  81.  * construct histogram of input data
  82.  */
  83. void
  84. construct_hist (int n, float *data, float *hist, double bw, double base)
  85. {
  86.   int i, ib;
  87.  
  88.   for (i = 0; i < n; i++) {
  89.     if (*(data + i) != 0) {
  90.       ib = 0;
  91.       while (*(data + i) >= base + ib * bw)
  92.         ib++;
  93.       *(hist + ib - 1) += 1;
  94.     }
  95.   }
  96. }
  97.  
  98.  
  99.  
  100. /*
  101.  * construct histogram of sorted(!) input data
  102.  */
  103. void
  104. construct_shist (int n, float *data, float *hist, double bw, double base)
  105. {
  106.   int i, ib;
  107.  
  108.  
  109.   ib = 0;
  110.   for (i = 0; i < n; i++) {
  111.     if (*(data + i) != 0) {
  112.       while (*(data + i) >= base + ib * bw)
  113.         ib++;
  114.       *(hist + ib - 1) += 1;
  115.     }
  116.   }
  117. }
  118.  
  119.  
  120.  
  121.  
  122. /*
  123.  * determine mean of area histogram in form of centroid:
  124.  * (ref: F. Family & P. Meakin, Phys. Rev. A40, 3836-3854 (1989))
  125.  */
  126. double
  127. find_hist_mean (int nb, float *hist, double bw)
  128. {
  129.   int ib;
  130.   double ssum, sssum, hi;
  131.   double mean;
  132.  
  133.   ssum = sssum = 0.0;
  134.   for (ib = 0; ib < nb; ib++) {
  135.     hi = (double) (*(hist + ib));
  136.     ssum += ib * hi;
  137.     sssum += ib * ib * hi;
  138.   }
  139.   if ((-1.0e-10 < ssum) && (ssum < 1.0e-10))
  140.     mean = -1.0;                /* error */
  141.   else
  142.     mean = (sssum / ssum) * bw;
  143.  
  144.   return (mean);
  145. }
  146.  
  147. /*
  148.  * prepare to evaluate histogram
  149.  */
  150. struct Histo *
  151. eval_hist (int BIN_METHOD, float *data, int nd, int SRT)
  152. {
  153. #ifdef HIST_DEBUG
  154.   int ib;
  155.   int id;
  156. #endif
  157.   int nb;
  158.   double bw, base, range;
  159.   struct Histo *h;
  160.  
  161.   switch (BIN_METHOD) {
  162.   case DEFAULT:
  163.     nb = nd + 1;
  164.     base = *(data + 0);
  165.     range = *(data + nd - 1) - *(data + 0);
  166.     if ((bw = range / (double) nb) < 1.0) {  /* for area data */
  167.       bw = 1.0;
  168.       nb = (int) range;
  169.     }
  170.     break;
  171.   case FIX_BINWD:
  172.     bw = (double) BIN_WIDTH;
  173.     base = *(data + 0);
  174.     range = *(data + nd - 1) - *(data + 0);
  175.     while ((nb = (int) (range / (int) bw)) > NBIN_LIM)
  176.       bw++;
  177.     break;
  178.   case FIX_RANGE:
  179.     nb = N_BINS;
  180.     base = (double) MIN_AREA;
  181.     range = (double) MAX_AREA - base;
  182.     if (MAX_AREA < *(data + nd - 1)) {
  183.       printf ("...WARNING: ");
  184.       printf ("max area exceeds currently set upper limit!\n");
  185.     }
  186.     bw = range / (double) nb;
  187.     break;
  188.   case FIX_INTERVAL:
  189.     nb = N_BINS;
  190.     base = 0.0;
  191.     range = 2.5;
  192.     if (MAX_AREA < *(data + nd - 1)) {
  193.       printf ("...WARNING: ");
  194.       printf ("max area exceeds currently set upper limit!\n");
  195.     }
  196.     bw = range / (double) (nb - 1);
  197.     break;
  198.   default:
  199.     printf ("\n...unknown binning method\n");
  200.     return ((struct Histo *) NULL);
  201.     break;
  202.   }
  203.  
  204.   if ((h = (struct Histo *) calloc (1, sizeof (struct Histo))) == NULL)
  205.       fail_alloc ("h", 1);
  206.  
  207.   if ((h->ph = (float *) calloc (nb, sizeof (float))) == NULL)
  208.       fail_alloc ("h->ph", 1);
  209.  
  210. #ifdef HIST_DEBUG
  211.   printf ("...evaluate histogram: nb = %d, bw = %f\n", nb, bw);
  212.   printf ("...input data\n");
  213.   for (id = 0; id < nd; id++) {
  214.     printf ("  %6.2f", *(data + id));
  215.     if ((id + 1) % 8 == 0)
  216.       printf ("\n");
  217.   }
  218.   printf ("\n");
  219. #endif
  220.  
  221.   if (SRT == 1) {               /* inp data sorted */
  222.     construct_shist (nd, data, h->ph, bw, base);
  223.  
  224.     h->amin = *(data + 0);      /* min, max inp vals */
  225.     h->amax = *(data + nd - 1);
  226.   }
  227.   else {
  228.     construct_hist (nd, data, h->ph, bw, base);
  229.  
  230.     h->amin = (float) -1.0;     /* min, max inp vals */
  231.     h->amax = (float) -1.0;
  232.   }
  233.  
  234.  
  235.  
  236. #ifdef HIST_DEBUG
  237.   printf ("...histogram\n");
  238.   for (ib = 0; ib < nb; ib++) {
  239.     printf (" %4.1f ", h->ph[ib]);
  240.     if ((ib + 1) % 10 == 0)
  241.       printf ("\n");
  242.   }
  243.   printf ("\n");
  244. #endif
  245.  
  246.  
  247.   h->nh = nd;
  248.   h->phd = &data[0];
  249.   h->mean = -1.0;
  250.   h->sdev = -1.0;
  251.   h->skew = -1.0;
  252.  
  253.   h->meth = BIN_METHOD;
  254.   h->nb = nb;
  255.   h->bw = (float) bw;
  256.   h->hmean = find_hist_mean (nb, h->ph, bw);
  257.  
  258.   return (h);
  259. }
  260.